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Abstract 

We are interested in simulating blood flow in arteries with a one di- 
mensional model. Thanks to recent developments in the analysis of hy- 
perbolic system of conservation laws (in the Saint- Venant/ shallow water 
equations context) we will perform a simple finite volume scheme. We 
focus on conservation properties of this scheme which were not previously 
considered. To emphasize the necessity of this scheme, we present how a 
too simple numerical scheme may induce spurious flows when the basic 
static shape of the radius changes. On contrary, the proposed scheme 
is "well-balanced": it preserves equilibria of Q = 0. Then examples of 
analytical or linearized solutions with and without viscous damping are 
presented to validate the calculations. The influence of abrupt change of 
basic radius is emphasized in the case of an aneurism. 

Keywords blood flow simulation; well-balanced scheme; finite volume scheme; 
hydrostatic reconstruction; man at eternal rest; semi-analytical solutions; shal- 
low water 

1 Introduction 

As quoted by Xiu and Sherwin [86] the one dimensional system of equations 
for blood flow in arteries was written long time ago by Leonhard Euler in 1775. 
Of course, he noticed (see Parker [64]) that it was far too difficult to solve. 
Since then, the pulsatile wave flow has been understood and is explained in 
classical books such as Lighthill [50] and Pedley [65]. More recently a huge lot 
of work has been done and a lot of methods of resolution have been proposed 
to handle this set of equations, among them [81], [88], [73], [62], [44] [42] [71], 
[86], [25], [83], [58], [27], [15], [1], [16], [84], [84, 55]... and we forgot a lot of 
actors. Even if recent progress of fluid structure interaction in 3D have been 
performed (Gerbeau et. al. [24], Van de Vosse et. al. [79], among others), those 
advances have not made the 1-D modelisation obsolete. On the contrary, it is 
needed for validation, for physical comprehension, for boundary conditions, fast 

* Corresponding author: DelestreOunice.fr, presently at: Laboratoire de Mathematiques 
J. A. Dieudonne (UMR CNRS 7351) - Polytech Nice-Sophia , Universite de Niee - Sophia 
Antipolis, Pare Valrose, 06108 Nice eedex 02, France 

tCNRS & UPMC Universite Paris 06, UMR 7190, 4 place Jussieu, Institut Jean Le Rond 
d'Alembert, Boite 162, F-75005 Paris, France 



or real time computations and it is always relevant for full complex networks 
of arteries ([88], [73], [69], [62, 63]) or veins ([28]) or in micro circulation in the 
brain [67]. 

In most of the previous references it is solved thanks to the method of char- 
acteristics, to finite differences or finite volumes (that we will discuss and im- 
prove in this paper), or finite elements methods (Galerkin discontinuous). Each 
method has its own advantages. 

In the meanwhile, the set of shallow water equations (though younger, as 
the equations have been written by Adhemar Barre de Saint- Venant in 1871 
[18]) have received a large audience as well. They are used for modelling the 
flow in rivers [31], [12] (and networks of rivers), in lakes, rain overland flow 
[20, 21], [75], in dams [77, 78] or the long waves in shallow water seas (tides in 
the Channel for example [70], or for Tsunami modelling [23], [68]). It has been 
observed that Saint- Venant equations with source terms (which are the shape of 
the topography, and the viscous terms) present numerical difhculties for steady 
states. The splitting introduced by the discretisation creates an extra unphysical 
flow driven by the change of the topography. A configuration with no flow, for 
example a lake, will not stay at rest: thus, the so called equilibrium of the "lake 
at rest" is not preserved. So new schemes, called "well-balanced" have been 
constructed to preserve this equilibrium. Bermudez and Vazquez [6, 5] have 
modified the Roe solver to preserve steady states. Other ways to adapt exact or 
approximate Riemann solver to non-homogeneous case have been proposed by 
LeVeque [47] and Jin [37]. Following the idea of the pioneer work of Greenberg 
and LeRoux [32], [30] and [29] have proposed schemes based on the solution 
of the Riemann problem associated with a larger system (a third equation on 
the bottom topography variable is added). A central scheme approach has 
been developped in [41]. In Perthame and Simeoni [66], the source terms are 
included in the kinetic formulation. In [3] and [8], a well-balanced hydrostatic 
reconstruction has been derived. It can be applied to any conservative finite 
volume scheme approximating the homogeneous Saint- Venant system. In [82], 
[17], [8] and [52, 51] the friction source term is treated under a well-balanced 
way. More recently, in [45], the system of shallow water with topography is 
rewritten under a homogeneous form. This enables them to get Lax-WendrofF 
and MacCormack well-balanced schemes. With well-balanced methods, spurious 
current does not appear. 

Thus we will follow those authors and construct by analogy a similar method 
but for flows in elastic tubes rather than for free surface flows. Even if for artery 
flows the wave behaviour of the solution is very important (forced regime with 
lot of reflexions), we have in mind complex networks. In this case we will 
have to deal from the arteries to the veins with eventually microcirculation 
in the brain. In the smaller and softer vessels, the ID approximation is still 
valid, the pulsatile behaviour is less important, and there are lot of changes of 
sections (across stenosis, aneurisms, valves). Those changes of sections have to 
be treated numerically with great care, which has not been the case up to now 
in the litterature. 

So here, we will follow the authors of the Saint- Venant community and con- 
struct by analogy a similar method but for artery flow. At first we derive the 
equations written in conservative form, looking in the literature, we found that 
it has been only written by [83] in a PhD thesis and by [15], [16] but not ex- 
ploited in the "well-balanced" point of view. Mostly ([86], [55]), equations are 
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in fact written in a non-conservative form. With non- conservative schemes, we 
have consistency defaults and shocks are not uniquely defined, see among others 
[36], [7] and [13]. Furthermore, the discharge Q is not conserved. Not a lot of 
attention has been attached to this problem up to now. But this non conser- 
vation is a real problem if we consider networks with tubes with non uniform 
section and ultimately quasi steady flow in small terminal vessels. 

Hence after the introduction of the set of equations in section 2 and deriving 
with this new point of view the conservative equations, we construct by analogy 
the numerical scheme both at first and second order in section 3. Then section 
4 is devoted to several test cases to validate the different parts of the equations. 
The most important case is the "man at eternal rest" by analogy to the "lake 
at rest" which would not be preserved by non-conservative schemes. The aim 
of this paper is to use this property to construct efficient methods of resolution. 



2.1 Derivation of the equations 

We first derive the model equations, they are of course now classical, as already 
mentioned. Starting from unsteady incompressible axi-symmetrical Navier- 
Stokcs equations, and doing a long wave approximation one obtains a first set 
of Reduced Navier-Stokes equations ([42], [43]). One has the longitudinal con- 
vective term, the longitudinal pressure gradient, and only the transverse viscous 
term, the longitudinal being negligible. The pressure is constant across the cross 
section and depends only of x the longitudinal variable. The incompressibility 
is preserved. 

Then, integrating over the section, one obtains two integrated equations 
(kind of Von Karman equation) . The first one is obtained from incompressibility 
and application of boundary condition and relates the variation of section to the 
variations of flow. The second one is obtained from the momentum. The final 
equations are not closed and one has to add some hypothesis on the shape of 
the velocity profile to obtain a closed system (some discussions of this are in 
[42]). We then obtain the following system with dimensions which is the 1-D 
model of flow: 



where A{x,t) is the cross-section area {A — ttR^, where R is the radius of the 
vessel), Q{x,t) = A{x,t)u{x,t) is the flow rate or the discharge, u{x,t) the 
mean flow velocity, and p the blood density (see figure 1). The first equation is 
without any approximation. An extra factor of value 4/3 may appear in front 
of Q'^/A if the chosen profile is a Poiseuille profile (with either large viscosity 
I' or low frequency uj / {2tt) so that a — Rq jv the Womersley number is 
small). This comes from the fact that this term is in fact an approximation 
of /q 2Tru{t,x,r)'^rdr (this u{t,x,r) being the actual longitudinal velocity, not 

to be confused with u{x,i) = 2'!Tu{t, x,r)rdr). The chosen unit value cor- 
responds to a flat profile. The viscous term has been modelised by a friction 
proportional to Q/A (see [87] [42]): C/ = Sni^ where v is the blood viscosity. 
Again, a friction term written —SnvQ/A is an approximation of the effective 
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term which is vRdu/ dr\fi (with here again u{t,x,r) the real longitudinal ve- 
locity). This specific choice corresponds to a Poiseuille flow (small a), and is 
therefore non consistent with the coefficient in the convective term. This is 
not a real problem, those coefficients are adjusted in the literature. Note that 
in Saint- Venant community, the skin friction modelling involves a square of Q, 
reminiscent to the turbulent nature of river flow. This is seldom used in blood- 
flow. To conclude on the mechanical modelling we will take the form of Eq. 1 
for granted and we will no more discuss it. 

Finally, one has to model the arterial wall, the pressure has to be linked with 
the displacement of the vessel. A simple elastic linear term in the variation of 
radius is taken: 

p — Po = k{R — Rq) or as weU p = po + k^^——7^—^, (2) 

Ao{x) is the cross section at rest {Aq = ttRq, where Rq is the radius of the 
vessel which may be variable in the case of aneurism, stenosis or taper) , and the 
stiffness k is taken as a constant. The external pressure po is supposed constant 
as well. Again, this relation which is here a simple elastic string model may 
be enhanced. Non linear terms, non local (tension) terms and/or dissipative 
unsteady terms may be introduced (viscoelasticity). We take here the most 
simple law, but we are aware of the possibilities of complexity. 

Note that the derivation of those equations may be done with only ID fluid 
mechanics arguments ([67] or [39]). One has to be very careful then to write 
precisely the action of the pressure on the lateral walls (for example Kundu [40] 
page 793 does a mistake in his text book). But, the real point of view to derive 
the equations is the one explained at the beginning of this section and involves 
the Reduced Navier-Stokes equations ([42]). 

We will now change a bit this form to write a real conservative form. Note 
that authors ([71], [55]) write system (1) as 

dtU + d^F{U) = S{U), 

they use the variables U = {A,u), so that they identify a kind of flux F{U) — 
{Au, u^/2+p/ p) and a source term S{U) ~ —Cfu/A. But u is not a conservative 
quantity, Q is, as we will see next. 

Of course, those Saint- Venant and blood flow equations are very similar to 
the ID Euler compressible equations that one can write in nozzles, or in acoustics 
(Lighthill [50]), in this case the density p is variable and linked to the pressure 
by the isentropic law p on p'^ . This relation is the counter part of the relation of 
pressure which are in Saint- Venant and respectively in the artery: hydrostatic 
balance p = pgrj (density, gravitational acceleration, level of the free surface) 
and respectively the elasticity response of the artery p = k{R — Rq) (stiffness of 
the wall, current radius of the vessel, reference radius at rest). 
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Figure 1: ID models: left the blood flow model {u is the mean velocity, R the 
radius, A = ttR^ the instantaneous area, Rq the radius at rest, Aq = ttRq the 
area with no flow), right the shallow water model (u is the mean velocity, h the 
water column height, z the topography). 



2.2 Conservative system 

So, to write the conservative system, the second equation of (f ) is rewritten and 
developed. Thus we get the following system of conservation laws: 



dtA + d^Q=^0 
dtQ + dx ' 



A 3p^/t: 

we can write (3) under a more compact form: 

dtU + d,F{U) = S{U), 



(4) 



with a new definition of U, F{U) and S{U). We identify U the vector of the 
conservative variables, F{U) the flux 



U = 



A 
Q 



Q 

F{U) = I ^^3/2 

A 3pv^ 



(5) 



and the source term takes into account the initial shape of the vessel Aq and 
the friction term: 



Py/TT A 



(6) 



This system is strictly hyperbolic when ^ > (which should be the case in 
arteries). Indeed, this means that we have 



1 



A 
Q 



= J{u).d^u, 



where the Jacobian matrix J{U) has two real eigenvalues (which is the definition 
of hyperbolicity) 



Ai = 



Q 

A 



IkVA A ^ 

= = u — c and A2 = — + 

2py/'K A 



I kVA 



= u + c. 



(7) 



We recognise c as the well-known Moens Korteweg wave propagation speed 
([65]): 

/ kVA 
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The system (3) without friction i.e. with Cf — admits special solutions which 
are steady state solutions, that means that we have 

dtA = dtu = dtQ = 0. (8) 

Thus the mass conservation equation gives the discharge conservation: dxQ{x,t) - 
0. The momentum conservation equation, under the steady state and frictionless 
hypothesis reduces to 



2 



0, (9) 



with b = fc/(/9Y^). So that system (3) integrates to the two constants 

,t) = Qo 

+ bVA - by/Ao = Cst 



Q{x,t) = Qo 

91.. n. .nr- • (10) 

2A^ 



The first is a constant flux, the second is the Bernoulli constant. 

We have of course made analogy between (3) and the one-dimensional Saint- 
Venant (or shallow water) system: 

dth + dxq^O 

where h{x,t) is the water height, u{x^t) the mean flow velocity, q{x,t) — 
h{x,t)u{x,t) the flow rate and z{x,t) the topography of the bottom and c is 
the friction coefficient. With /3 = 2 we recover Chezy's or Darcy-Weisbach's law 
(depending on the way c is written) and with /3 = 7/3, we get the Manning- 
Strickler's friction law. 

This system admits steady state frictionless solutions too: 

dth dtU = dtq = 0. (12) 

Thus we get the conservation of discharge and the Bernoulli's relation 

9 = go 

2gh^ 

which is the exact analogous of (10). In the literature we can find numerical 
schemes preserving steady states under the form (13), but they are complex to 
handle (see [14], [60], [76] and [9]). Most of the time, more simple steady states 
are considered 

q = u = 

dx[9^^+ghdxZ = Q ' (1^) 

which correspond from a mechanical point of view to the modelling of a " lake 
at rest" . We have a balance called hydrostatic balance between the hydrostatic 
pressure and the gravitational acceleration down an inclined bottom z. The 
second term of (14) reduces to 

H = h + z = Cst, (15) 
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where H is the water level. Since the work of [32], schemes preserving at least 
(14) are called well-balanced schemes. This allows to have an efficient treatment 
of the source term. Thus, in practice, schemes preserving (14) give good results 
even in unsteady cases. The analogous of the "lake at rest" for the system 
(3) without friction, can be called the "man at eternal rest" or "dead man 
equilibrium" , it writes 

In this case we have vC4 = ^/ Aq. We will now use this property to construct 
the schemes. 



t 

n+1 
t" 








X 




















\ 




\ 









1/2 3^1 


1/2 ^i^l 



Figure 2: The time and space stencil. 



3 The numerical method 

In this section, we describe the scheme for system (1) at first and second order, 
with the space and time discretisation illustrated in Figure 2. 

3.1 Convective step 

For the homogeneous system i.e. system (4) without source term: S'(C/) = 0, a 
first-order conservative finite volume scheme writes simply (see figure 2): 

At Ax ' 

where U" is an approximation of U 



Ax 



-1/2 



n refers to time t„ with t„+i —tn — At and i to the cell Ci ~ {xi-l/2^ a^i+1/2) — 
{xi_i/2^ Xi_i/2-\' Ax). The two points numerical flux _F'j_|_]^/2 is an approximation 
of the flux function (5) at the cell interface z -I- 1/2 

F,+^,2=F{U,,U^+^). (18) 

This flux function will be detailed in section 3.3. 
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3.2 Source term treatment 
3.2.1 "Topography" treatment 

kA 

The source term — i=dx\/Ao is the analogous of the topography source term 

{—ghdxz) in the shallow water system (11). Looking at (3), it might be treated 
explicitly which writes 

ur' = (^:r+i/2 - ^11/2) + ^tsiun m 

where the source is written simply by evaluating the derivative of Aq: 

I 

S{U^)=\ kAf ^o.,+i/2-^»-i/2 

V 2pv^v^ Ace 

with the definition 



A, 



Oi + l/2 



2 

We will illustrate the fault of this naive method in sections 4.4 and 4.5. We 
will prefer a well-balanced method inspired from the hydrostatic reconstruction 
(see [3] and [8]). As for the hydrostatic reconstruction we do not consider the 
friction term and we are looking for a scheme which preserves steady states at 
rest (16), the second equation in (16) means that we have 

y/A - yj% = Cst, (20) 

as for the hydrostatic reconstruction we will use locally this relation to perform 
a reconstruction of the variable A. On each part of the interface, we have the 
following reconstructed variables 



y/A+i/2L = max(V37 + min(AVA^i+i/2,0),0) 

Ui+l/2L — {^i+l/2L^ ^i+l/ 2LUiy 

\/A+i/2B. = niax(y^A,+i - max(AV^j+i/2,0),0) ' 

, Ui+1/2B. = (^j + l/2_Ri ^i+l/2_RWi+l)* 



(21) 



where = V^o^+i - V^t- 

For consistency the scheme (17) has to be slightly modified under the form 



u^' - c/r - |^(^r+i/2L - F-_,/,^), (22) 

(23) 



where 



with 



Fi+1/2L - ^1+1/2 + Si+1/2L 



K"-1/2R - Fi-1/2 + ^i-l/2R 



Si+l/2L — 
Si-1/2R = 





V{Al^)-V{A^^,,^^) 


V{A-)-V{AU/2b) 



(24) 
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and 

V{A) = —^^3/2. (25) 

The numerical flux ^7+i/2 '^'^^ calculated with the reconstructed variables 
(21): 

Having treated the term linked to the variation of radius we now turn to the 
friction term. 

3.2.2 Friction treatment 

Concerning the friction treatment two methods are presented here, namely: a 
semi-implicit treatment ([26], [11], and [49]) which is classical in shallow water 
simulations and the apparent topography method (see [8], [54, 53]) which is a 
well-balanced method. 

Semi-implicit treatment We use (22) as a prediction without friction, i.e.: 

and a friction correction (via a semi-implicit treatment) is applied on the pre- 
dicted variables {U*): 



--Cfur\ (27) 



At 

thus we get u"^^ and we have A"^^ = ^i- H preserves zero velocity. 

Apparent topography The apparent topography method (see [8]) consists 
in writing the second equation of the system under the following form 



rr)2 



dtQ + 

with 



k 



A 3p^ 
dx\jAoapp = dx (y^) + bj 



-—j^dx^ Aoapp, (28) 



where _ 



d-rh ■■ 



k A^' 

The hydrostatic reconstruction (21) is performed with the corrected term y/Aoap-^ 
3.3 Numerical flux 

Several numerical fluxes might be used, some of them are defined in the follow- 
ing. 
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Rusanov flux Following [8], the Rusanov flux writes 

F{Ul)+F{U„) Un - Ul 
^(Ul, Uji) = c , 

with 

c= sup ( sup \\j{U)\), 
u=Ul,Ur ie{i,2} 

where Ai(C/) and X2{U) are the eigenvalues of the system. 

HLL flux Following [8] , the HLL flux (Harten, Lax and van Leer [35] ) writes 

FiUL) ifO<ci 

^[Ul,Ur) = < \ (Ur-Ul) ifci<0<C2 , 

F{Ur) if C2 < 

(29) 

with 

•^1 ^ ( inf Aj(C/))andc2= sup ( sup Xj{U)), 

U=UL,Ua je{l,2} U=Ul.Ur je{l,2} 

where Ai(C/) and X2{U) are the eigenvalues of the system. 

VFRoe-ncv flux Adapting the main line of [56] and [57], we get the following 
VFRoe-ncv flux for system (3) 

F([/L)ifO<^Ai(C/) 
J'{Ul,Ur)^{ F{U,)ii\i{U)<0<X2iU) , 
F{Ur) if A2(C/) <0 



with the mean state 



2 (CL + cr) 
W =[ ] = \ ul + ur 



Ac 



and 



A, 
A^ 



where et are the Roe mean states, we get them thanks the Riemann 
invariant 

thus we have 

( = u - 2{cr - cl) 
\ c^=c - ^{ur-ul) 

When Xi{Ul) < < Xi{Ur) or X2{Ul) < < X2{Ur), we can get non entropic 
solutions. We perform an entropy correction thanks to the Rusanov flux. 
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Kinetic flux Adapting the main line of [2] and [4] , we get the following kinetic 
flux for system (3) 

^(C/l, Ur) = F+iUL) + F^{Ur), 



with 



F+iUL) 



2V3TE 



Ah 



Ah-' 



F^Ur) 



Ar 



TO-)- 



where 



ith 



?>P\/tt 



A'h = max(0,UL 
M_ = max(0, 
TO-i- = min(0, Ur - 
TO_ — min(0, Ur - 

''Al and To — 



CFL condition We have to impose a CFL (Courant, Friedrichs, Levy) con- 
dition on the timestep to prevent blow up of the numerical values and to ensure 
the positivity of A. This classical stability condition writes 



At < ncFL- 



min(Aa:i) 



where Ci — a/ k\[A^il(2p^^ and ticfl — 1 (resp. ncpL — 0.5) for the first 
order scheme (resp. for the second order scheme see 3.5). 

We have to notice that the kinetic flux needs a particular CFL condition 
(see [2, p.24]) 

min(Aa;i) 

At < ncFL ' I r^. - 

max(|ui| + ^jAli) 

i 

The Rusanov is the most diffusive but the most simple to implement, the 
kinetic one is slightly less diffusive but more cpu time consuming. VFRoc-ncv 
and HLL are comparable, but the first one is a little more time consuming due 
to entropy correction. All these fluxes are compared in [19] in Saint- Venant 
framework. 



3.4 Boundary conditions 
3.4.1 Characteristic system 

Of course boundary conditions are very important for artery flow. We will 
not too much insist on them, and we will not present for example the classical 
Windkessel model or variations. Nevertheless, adapting [11] to blood flow, we 
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can write the homogeneous form of (4) under the following non conservative 
form 

dtA + udxA + AdxU = 

dtu H i=dxA + udxU — ' ^"^^^ 

2p^^/A 

Thanks to the Moens Korteweg velocity rewritten as A~ c*7r(2p/fc)^, we get 

dt [4c - it] + (u - c)5^ [4c - u] = 
dt [4c + m] + (u + c)dx [4c + w] = ' ^ ' 

d(Ac — u) d(A:C + u) 

thus we have = (respectively — 0) or the Riemann in- 

dt dt 
variant 4c — u = Cst (resp. 4c + u = Cst) along the characteristic curve C_ 

dx dx 
(resp. C-i_) defined by the equation — u — c (resp. — = u + c). 

dt dt 
The boundary conditions will be prescribed thanks to the method of char- 
acteristics. We set Ubound = U{0) or U{L) and Unear = U{Ax) or U{L — Ax) 
depending on the considered boundary (x — or x — L). 

3.4.2 Subcritical flow 

We write here boundary condition for subcritical flow: at the boundary the flow 
is subcritical if we have [uboimrfl < Cbound or equivalently 

{Uhound Abound) {ubound ^~ Abound) ^ 0- (33) 

Of course it seems to be always the case in blood flows, where |Q/(Ac)| is 
less than 10 % in physicological cases. This concept is more relevant to Saint- 
Venant equation, where supercritical flow are easily obtained. Two cases are 
possible: we are given either the cross section A (the pressure thanks to (2)) or 
the discharge Q. 

• Given cross section Abound = Cst: 

At a; = the Riemann invariant is constant along the ingoing characteristic 
u ~ c, thus we have 

Abound = Cst ^ ^ (34^ 

Ubound U^iear ^ linear ^bound\ 

At x = L the Riemann invariant is constant along the outgoing charac- 
teristic u + c, thus we have 

Abound — Cst (35) 
Ubound — ^near ~t- 4 [Cnear Abound] 

If (33) is not verified by the couple {Abound, Ubound), the flow is in fact 
supercritical. 

• Given discharge Qbound = Cst: 
Depending on the boundary considered we have either at a; = 

— Qbound ~t- {Ujif^ar ^^near) Abound ~t- ^Cbound Abound 7 (36) 
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or at X = L 



We recall that ctound depends on Abound, we solve (36) or (37) to get 

A-bou7id- 

The source terms might be considered negligible in the characteristic method. 
3.4.3 Supercritical flow 

As previously mentioned, this is not really a relevant case for flows in arteries. 
But, in veins, or in collapsible tubes, it may be relevant (see Pedley [65]). A 
supercritical inflow corresponds to both Abound — Csti and Qhound — Cst2 that 
have to be imposed. A supercritical outflow is such that we have Abound = Anear 

and Q bound Qnear- 

Again, source terms are neglected in the characteristic method. 

In order to impose the discharge Qhound = AhoundUbound, we can use an other 
method: the use of the first component of the numerical fiux is an accurate way 
to proceed, even if J^i{Ul,Ur) — Qbound{t) has no unique solution (and the 
complexity of the problem depends on the numerical flux). 

3.5 Second order scheme 

In order to get a second order scheme, we use the following algorithm: 

Step 1 In order to get second order reconstructed variables (C/,±, VAo»±) we 
perform a linear reconstruction on variables u. A, \fA — -s/A] then we get the 
reconstruction on \fAQ. 

We consider the scalar function s € M, its linear reconstruction is defined by 

Aa; Ax 

Si-l/2+ = Si 2~^*' = Si + -^DSi 

where D is one of the following reconstruction operator Dmusch D^no and 
Denom- To get the reconstruction operator D, we introduce the minmod slope 
limiter 

{min(x, y) H x,y > 
ma.x{x,y) iix,y<0 . 
else 

Some other slope limiters might be found in [48, p. 111-112] (the Monotonized 
Central difference limiter and the Superbee limiter). 



MUSCL With the operator 

. f Si ~ Si_l Si+l — Si \ 

DmusciSt = mmmod I — — — , — — — 1 , 

we get the MUSCL linear reconstruction (Monotonic Upwind Scheme for Con- 
servation Law [80]). We can find a less diffusive but more complicated form of 
the MUSCL reconstruction in [22]. 
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ENO With the operator 

DenoSi = minmod I — h Oeno-^D Si_i/2, h Oeno-^D Si+i/2 

where 



^ ^«+i/2 = mmmod ( , 



and 



we perform the ENO hnear reconstruction (Essentially Non Oscillatory [34, 33], 
[72]). This reconstruction is more accurate than the MUSCL reconstruction but 
less stable. In order to get a more stable reconstruction we can take 9eno in [0, 1] 
(it is recommended to take 9eno = 0.25 [10, p.236]). 

modified ENO With the operator 

DenomSi = minmod {DenoSi, '20eno7nDrnusclSt) , with 0enom G [0, 1], 

we get the modified ENO linear reconstruction [8, p. 55-56] and [10, p. 235-237]. 
The previous linear reconstructions applied to A writes 

A,_i/2+ = A,~ ^DA„ A,+y-2- = A, + ^DA„ 
thus we have the following conservation relation 

2 ~ 

In order to get the same kind of relation on the discharges, the following modified 
reconstruction is performed on the velocity variables u 

Aj+i/2_ Aa; A,_i/2+ Ax 
Ui-i/2+ =Ui -J ^Du,, u^+i/2- =Ui + — — ^Du„ 

this allows to get the conservation relation 

^i-l/2+'"j-l/2+ + Ai+l/2-'"i+l/2- _ . 

2 

We perform on ^ = \/A — y/Aa the same linear reconstruction as on A, we get 
and at last we have 



V ^0,;+i/2± - Y "^«+l/2± ^ *i+l/2±- 

Step 2 From variables U,- and U,+ we get the following reconstructed vari- 
ables 

y^A^^i^jJ^ = max(^'j+i/2_ + mm{y^,+i/2-,VA)i+i/2+),0) 

Ui +1/2L = (Ai+1/2L, A'i+l/2L"'i+l/2-)* /oo^ 

y/A+i/2R = max(5'i+i/2+ + min(VA)»+i/2-, VA)^+i/2+), 0) 

Ui+1/2R — iAi+i/2R, Ai^if2RUi+i/2+y 
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step 3 In order to get a consistent well-balanced scheme, it is necessary to 
add a centered source term Fci. Scheme (17) becomes 



where 



with 



FT+1/2L - ^1+1/2 + -5*4+1/2- 



Si+l/2- 





^(^m/2J-^(^r+i/2L) 



^'-^/^+~(^(^r-i/2+)°^(^r-i/2i.,) 



and 







with Ay/Aoi = •\/A)i+i/2_ — •\/A^i_i/2+- The numerical flux -F'j"j_i/2 calculated 
as at first order (3.2.1), i.e. with the variables obtained with the hydrostatic 
reconstruction (38). 

Thus we get a second order scheme in space which writes 

f/"+i = U" + Ai$(C/") with U = ([/i)iez 

and 

mn = {Fr+i/2L - Fr-i/2R - Fdi) . m 

step 4 The second order in time is recovered by a second order TVD (Total 
Variation Diminishing) Runge-Kutta (see [72]), namely the Heun method 

fjn+l ^ jjn ^ At$(C/"), 

jjn+2 ^ jjn+1 ^ A<$(C/"+l), /^QN 

2 

where $ is defined by (39). This method is a kind of predictor-corrector method. 

A more complete modelling of the arterial flow will imply other source terms 
(diffusion, tension,...). To observe those effects one has to use a higher order 
scheme, as performed for shallow water in [59] and [57]. 



4 Validation 

4.1 About the chosen examples 

The chosen examples are more for sake of illustration of the previous scheme 
rather than for validation on real clinical datas. So, the extremities of the 
arteries will be non reflecting and the numerical values are more indicative rather 



15 



than cHnically relevant. Those examples will show that the scheme behaves well 
and that some validations from linearised or exact solutions may be reobtained. 
More specifically we insist on the cases with large change of initial section like 
sudden constriction and sudden expansion. Even we present here an aneurism, 
flow in a stenosis could also be treated. 

4.2 The ideal tourniquet 

In the Saint- Venant community, the dam break problem is a classical test case 
[74] (in compressible gas dynamics, it is referred as the Sod tube, LeVeque [46] 
or Lighthill [50] or [61]). We have an analytical solution of this problem thanks 
to the characteristic method. So, we consider the analogous of this problem 
in blood flow: a tourniquet is applied and we remove it instantaneously. Of 
course, this is done in supposing that the pulsatile effects are neglected, so it 
is more a test case rather than a real tourniquet. Any way, with the previous 
hypothesis, we get a Riemann problem. This test allows us to study the ability 
of the numerical scheme to give the front propagation properly, we notice that 
the non linear terms are tested (but neither the viscous ones nor the change 
of basic radius). As we consider an artery with a constant radius at rest and 




Figure 3: (a) The ideal tourniquet radius initial conditions: a discontinuous 
initial radius (or pression) is imposed with no velocity (the blood is blocked), 
(b) analytical solution of the radius at time t > 0: a shock wave is moving 
downstream (to the right) and an expansion wave upstream (to the left). 

without friction, we consider the system (3) without source term. L is the length 
of the arteria with x E [—L/2, L/2]. The initial conditions are 

with Al > All and Q(x,0) = ^(a;, 0) — (see figure 3). 

As for the dam break, with the characteristic method we get an analytical 
solution. At time t > 0, from the upstream to the downstream of the artery, 
we have four different states (as illustrated for the radius R = \/ A/n in Figure 
3b): 
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1. Upstream point A located at XA{t) = —CLt, the state is the same as the 
initial one: 

u{x,t)=0 • ^ ' 

2. Between point A and point B located at XB{t) = {um — cmY = (4ci — 
5cA/)i (the subscript stands for the intermediate state in the following 
zone), we have 

Ax 4 

(43) 

, , 1x4 
c{x,t) = + -Cl 

5 t 5 

3. Between point B and C located at xc{t) = st, with the Rankinc Hugoniot 
relation s — AmUm /{Am — Aj^). We have an intermediate state 

A{x,t)^AM ,44. 

U{X, t) = UM ' 

where the variables s, Am, "a/ and Qm — umAm are obtained thanks to 
the following system 

UM + 4CM =UL+ AcL 

Qr - Qm = s{Ar - Am) 
As Ufj = = Qr = Ql = 0, this system reduces to 

UM + ^CM = 4CL 

Qm = s[Am - Ar) 

Qm I ^ A 3/2] k A 3/2 _ 
^ Am ~ 7, — i=Aii — sQm 



Am Spy/TT J Spy/TT 

This system is solved iteratively. 

4. Downstream point C, the state is the same as the initial one 

A{x,t) = Ab 
u{x, t) = Q 



(45) 



The simulations were performed with the first order scheme with each of 
the fluxes, J = 100 cells and a fix CFL of 1 (0.5 for the kinetic flux) and 
we have used the following numerical values C/ = 0, fc = 1.0 lO^Pa/m, p = 
1060m3, Rl = 5. lO-^m, = 4. IQ-^m, L = 0.08m, T^nd = 0.005s, cr = 
yJkRR/{2p) ~ 4.34m/s, cl = y/kR^{2p) ~ 4.86m/s. An initial flow at rest: 
Q{x,t = 0) = OmVs and 



A{x,t = 0) 



ttRl\ if X G [-0,04 : 0] 
ttRr^, if a; e]0 : 0,04] 



As illustrated in Figure 4, we notice that the Rusanov and Kinetic fluxes are 
more diffusive than the others. With mesh refinement and/or scheme order in- 
creasing, we should improve the results both at the shock level and the expansion 
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Figure 4: Snapshots of the ideal tourniquet - comparison between analytical 
and numerical results at t = Tend- (a) radius and (b) discharge. 



wave. As noticed in [19] for the shallow water equations, the VFRoe-ncv needs 
more cpu time than the HLL flux. So in what follows we will focus on the HLL 
flux. This is a rather severe test case as in practice, |Q/(Ac)| is less than 10 % 
in physicological cases. 

4.3 Wave equation 




Figure 5: The propagation of an initial (non physiological) pulse given at t = 
in space and time: (a) radius and (b) velocity: the initial perturbation splits 
into a forward and backward travelling wave (colors online). 

So, as in practice, in arteries, the mean velocity is far smaller than the 
Moens Korteweg celerity |(5/(v4c)| ^ 1, here we will study the ability to catch 
the relevant small linearised wave propagation speeds. We look at the evolution 
of a given perturbation of radius with no velocity perturbation and see how it 
evolves. Thus we consider the system without friction, with a constant section at 
rest, a steady state (Rq, uq = 0). Looking at perturbation at a small level £ <C 1, 
we have {R,u) = {Rq + eRi^eui) where subscript 1 refers to the perturbation. 
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Substitution in (3) and neglecting the small terms, it becomes the d'Alembert 
equation for the perturbations ui and 

dUu,,R,)~co^dUu,,R,)=0, (46) 

the same is valid for variables Qi = AqUi the perturbation of flux and pi — kRi 
the perturbation of pressure with the Moens Korteweg wave velocity 



Co 




(47) 

[x) we obtain so analytical 



With initial conditions ui{x,Q) = 0, Ri{x,{)) - 
solution: 

' R[x, i) = i?o + I [0 (x - Cot) + (t){x + Cot)] 

u{x,t) = [-(j} {x - Cot) + (j) {x + Cot)] 

I Ho 

The following numerical values have been used for the figure 5: J = 200, Cf — 0, 
k = 1.0 10*^Pa/m, p = 1060m3, Rq = 4. lO^^m, L = 0.16m, T^nd = 0.008s, 
Co — \/kRo/{2p) ~ 13.7m/s. As initial conditions, we take a fluid at rest 
Q(x, 0) = Om'^/s with an initial deformation of the radius 



A{x,0) = 



7ri^o^ if x G [0 : 40L/100] U [60L/100 : L] 



■kRo 



e sm 



x-40L/100\ 
20L/100 J 



, if X e]40L/100 : 60L/100[ ' 



with e = 5.10"'^. As illustrated, in Figure 6a, we get two waves, propagating on 
the right (respectively left) with a positive (resp. negative) velocity. The two 
waves are represented at several moments. We notice no numerical diffusion. 
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y=-0.78x-17.03 






Order 2 






y=-1.15x-15.42 
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log(J) 
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(a) 



(b) 



Figure 6: (a) Radius as function of the position at 3 time steps: t — 0, t — 
Tend/^, t = Tend/ '2' and t = ST^nd/^, Comparison between the first order scheme 
with HLL flux and exact solution and (b) the error made on the area calculation 
at t = 0.004s. 
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4.4 "The man at eternal rest" 



The previous test case did not involve drastic changes in the basic radius. We 
will show in this section that a not adapted source term treatment in (3) may 
give non-physical velocity. In this test case, we consider a configuration with 
no flow and with a change of radius Rq{x), this is the case for a dead man with 
an aneurism. Thus the section of the artery is not constant and the velocity is 
u{x,t) = m/s. We use the following numerical values. J = 50 cells, C/ ~ 0, 
k — 1.0 lO^Pa/m, p — IGGOm"^, L = 0.14m, Tend = 5s. As initial conditions, we 
take a fluid at rest Q{x, 0) = Om-^/s and 



R{x,0) = Ro{x) 



i?o if X e [0 : xi] U [x4 : L] 

AR r 

Ro H — ^ sm 



^X2 - Xi 

Ro + AR if X e [x2 : X3] 
P ^ Ai? / f X-X3 

Ro + cos 



-7r-7r/2 +1 



, if X e]xi : X2[ 



if X g]x3 : X4[ 



with i?o = 4.0 IQ-^m, AR = 1.0 10"^ 



m, xi 



= 1.0 lO-^m, X2 = 3.05 lO-^m, 



X3 = 4.95 lO^^m and X4 = 7.0 lO^^m (figure 7a). 

As illustrated in Figure 7b, with an explicit treatment of the source term in 
Aq^ we get non zero velocities at the arteria cross section variation level. The 
"man at eternal rest" is not preserved, spurious flow appear (avoiding artifacts 
such as [38]). As expected, the hydrostatic reconstruction (21) preserves exactly 
the steady state at rest. 




Figure 7: The "dead man case": (a) the radius of the arteria, (b) comparison 
of the velocity at time t = 5s between an explicit treatment of the source term 
and the hydrostatic reconstruction. The spurious flow effect is clearly visible if 
an inappropriate scheme is chosen. 



4.5 Wave reflection-transmission in an aneurism 
4.5.1 Analytical linear transmission 

We observed that waves propagate fairly well in a straight tube, we now observe 
the reflexion and transmission trough a sudden change of section (from to 
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Ar) in an elastic tube. The d'Alembert equation (46) admits harmonic waves 
^t{ujt-x/co) ^ so that a plane wave is (symbol 1 represents a linear perturbation 
as in section 4.3): 

= Yopi, with pi oc 7ee(e*("*"^/"«') 

and Yq = Ao/{pco) is called the characteristic admittance. With this definition, 
we can look at a change of section from a radius i?^ to a radius Rji: an imposed 
right traveling plane wave g^i^^-^/'^i') moving at celerity cl in the left media 
(subscript L) will experience a transmission in the right media (subscript R) 
(and will move at celerity cn) and a reflexion (and move at celerity — c^). The 
coefficients of transmission and reflexion depend on the two admittances of the 
two media L and R (see Lighthill or Pedley): 

Tr = — and Re = — 

Yl + Yr Yl + Yr 

As the admittance does not depend on the frenquency w, any signal (by Fourier 
decomposition) will be transmitted and reflected with those values. 

Here we will study the reflection and the transmission of a small wave in an 
aneurism. 

4.5.2 Propagation of a pulse to and from an expansion 




Figure 8: To an expansion, the propagation of an initial pulse in space and time: 
(a) radius and (b) velocity (colors online). 



First we test the case of a pulse in a section Rr passing trough an expansion: 
Al > Ar, taking the following numerical values: J = 1500 cells, Cf — 0, 

AR = 



lOeOm'^ 



= 5. lO-^m, Rr = 4. 

-3c 



10 



-3, 



k = 1.0 lO^Pa/m, p 

1.0 10-^m, L = 0.16m, Tend = 8.0 lO'^s, cl = ^/kR^J(2p) ~ 15.36m/s and 
cr = yJkRR/{2p) ~ 13.74m/s. We take a decreasing shape on a rather small 
scale: 



Ro{x) 




X — Xi 
X2- Xi 



if X e [0 : xi] 
if X e]xi : X2] 
else 



(48) 
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eR(,Tr/2 












ERoRe/2 
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x[m] 



Figure 9: Radius as function of the position at 3 time steps: t = 0, t ^ Tend/4:, 
t = Tend/'^ and t = 3Tend/4:. The dotted hnes represent the level of the predicted 
reflexion (Re) and transmission (Tr) coefficients. 



with xi — 19L/40 and X2 = L/2. As initial conditions, we consider a fluid at 
rest Q{x, 0) — Om'^/s and the following perturbation of radius: 



R{x,0) = 




. /too 

e sm TT 

\20L 



65L 

Too 



if x e [65L/100 : 85L/f00] 
else 





Figure f 0: From an expansion, the propagation of an initial pulse in space and 
time: (a) radius and (b) velocity (colors online). 

Results showing the wave propagation and expansion are in Figure 8. In 
Figure 9 we have the amplitude of the pulse before and after the change of 
section, the level of the analytical prediction of T and R is plotted as well 
showing that the levels are preserved. 

The same is done for a pulse propagating from an expansion. So, just the 
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Figure 11: Radius as function of the position at 3 time steps: t — 0, t = Tend/4:, 
t = Tf-nd/'i and t = 3renci/4. The dotted lines represent the level of the predicted 
reflexion {Re) and transmission (Tr) coefficients. 



radius is changed: 



i?(a;,0) = 



Ro{x) 



. /lOO / 



15L 

Too 



if x e [15L/100 : 35L/100] 
else 



with e = 5.0 10"'^. Similar results showing the wave propagation are in Figure 
10a and 10b. In Figure 11 the amplitude of the pulse corresponds to the level 
predicted by the analytical linearised solution. 



4.5.3 The aneurism 
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(a) 



(b) 



Figure 12: The propagation in an aneurysm of an initial pulse in space and time: 
(a) radius and (b) velocity. We note the expected reflexions at the position where 
the vessel changes his shape (colors online). 



In this part we put an expansion and one constriction modelling an aneurism. 
An initial pulse at the center C of this aneurism propagates in the right and left 
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Figure 13: A scheme which does not preserve the "man at eternal rest" generates 
spurious waves were the radius changes (colors online). 



direction to the boundaries B. One sees in Figure 12 the reflexions and trans- 
missions at the position of the variation of radius. These plots have been done 
with the following numerical values: J = 1500 cells, Cf — 0, k ~ 1.0 10*Pa/m, 
p = 1060m3, Rg = 4. lO-^m, Rc = 5. lO^^j^, AR = 1.0 lO^^^^, L = 0.16m, 
T^d = 8.0 10-3s, CB = ^kRB/{2p) ~ 13.74m/s and cc = ^kRc/{2p) ~ 
15.36m/s. With the given shape 



Rb 
Rb 
Rb 
Rb 



1 



Ai? 

~Y 

AR = Rc 
AR ^ 



cos 



X 



Xi 



cos 



X2 - Xi 
X3 



X 



Xi - X3 



iix e [0 : xi] \J[xi : L] 
if X e]xi : X2] 

if a:: e]a;2 : 2:3] ' ^"^^^ 

if X e]x3 : X4[ 



with xi — 9L/40, X2 = L/4:, X3 = 3L/4 and X4 — 31L/40. As initial conditions, 
we consider a fluid at rest Q{x,0) = Om^^/s and the initial perturbation: 



i?(x,0) 



Roix) 
with e = 5.0 lO"^^ 



1 



- e sm 



100 

loZ' 



45L 



100) 



if .X e [45L/100 : 55i/100] 
else 



4.5.4 Non adapted treatment of the source terms 

To insist on the adapted treatment of the source terms, we present in this sub 
section what happens if a naive scheme is taken like in subsection 3.2.1. We 
clearly see in Figure 13 that at the change of radius reflected and spurious 
wave appear. The previous constant velocities are present like in Figure 13b. 
Furthermore, the initial data give traveling waves coming from the position of 
the change of section. This clearly shows the influence of the source terms in 
the discretization. 
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4.6 Wave "damping" 



In this last test case, we look at the viscous damping term in the linearized 
momentum equation. This is the analogous of the Womersley [85] problem, 
we consider a periodic signal at the inflow which is a perturbation of a steady 
state (i?o = Cst,UQ = 0) with a constant section at rest. We consider the set 
of equations 3 with the friction term under the nonconservative form with the 
variables (i?, u) , this system writes 



dtR + ud^R ^ 



R 



drU — 



dtu + udxU H — d^R 
P 



Cf u 



We take i? = i?o + sRi and u = + eui, where mi) is the perturbation of 
the steady state. We have after linearisation of the equations: 



2dtRi + Rod^ui 
k 

dtui + -dxRi = 
P 





^ Ro^ 



which may be combined in two wave equations with source term 

CfdtiuuRi) 



dt{ui,Ri)-co'di{ui,Ri) 



Rn 



(50) 



Looking for progressive waves {i.e. under the form e'*^"* ^^''), we obtain the 
dispersion relation 



ujC 



f 



Cq^ TtRqCo^ ' 

with uj = 2-K /Tpuisei Tpuise the time length of a pulse and K = k,. 
wave vector not to be confused with the stiffness, 



(51) 

ikj, the 



ki = 



w* 


( 




Ui?o'co2 




( 


Co* 


Ui?o'co2 



-I 1/4 



COS 



arctan — 



Cf 



ttRo'^u} 



1/4 



sm 



arctan — 



C 



f 



ttRo uj 



(52) 



For the numerical applications, we impose the incoming discharge 

Qbit) = Qamp sin(wt)m^/s, 

with Qarnp the amplitude of the inflow discharge. We get a damping wave in 
the domain 

^^^'^^ I Qamp sin (cji — fcr2^) e'^*'^ , '\ikrX<ujt ' ^^'^^ 

where Qbit) is the discharge imposed at x = 0. 

The following numerical values allow to plot figure 14a to 14d. For Cf, we 
took 4 diflFerent values {Cf = 0, C/ = 0.000022, Cf = 0.000202 and Cf = 
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Figure 14: The damping of a discharge wave (a.) Cf = 0, a — oo, (b) Cf — 
0.000022, a = 15.15, (c) Cf = 0.000202, a = 5 and (d) Cf = 0.005053, a ^ 1. 
The friction has been treated with either the semi-imphcit method (SI) or the 
apparent topography method (AT) and J = 800. 
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0.005053, corresponding to Womersley parameters a = oo, a ~ 15.15, a = 5 
and a = 1), J = 50, 100, 200, 400, 800 cells, k = l.lO^Pa/m, p = lOGOm^, 
Rq — 4.10^'^ni, L — 3ni, Tend — 25. s. As initial conditions, we take a fluid at 
rest Q{x,0) = Om'^/s and as input boundary condition 

Qb{t) = Qamp sin(w<)m^/s. 

with Lj = 2Tr/Tpuise — 27r/0.5s and Qamp = 3.45.10^^ni'^/s. As the flow is 
subcritical,the discharge is imposed at the outflow boundary thanks to (53) 
with X — L. 

In this part, we insist on the comparison between first order and second 
order. So, we compare first order (HHLl) and second order scheme with both 
MUSCL (HLL MUSCL) and ENO (HLL ENO) linear reconstruction. When 
C/ is not null, semi-implicit (SI) and apparent topography (AT) are compared. 
We should remark that as the friction increases the structure of the system 
(1) changes and goes from a transport/ wave behavior to a diffusive behavior. 
The diffusive behavior comes from the fact that as C/ increases, finally the 
momentum equation contains at leading order only the friction term and the 
elastic one: CfQi/Ao ^ —{AQk/p)dxR and so, using the conservation of mass: 

dR d^R . , ^ knRl 
- = 7.^w.thi..^, 

which is consistent with (51) at small w. However solving a "heat equation" with 
a "wave technique" seems to be adapted. Indeed, we notice a fair agreement 
between the numerics and the analytical solution in Figures 14a to 14d. "Errors 
made" between numerical results and the exact solution of the linearized system 
are given for information with CPU time in tables 1 and 2. 



4.7 Discussion on numerical precision 

Indeed, the notion of scheme order fully depends on the regularity of the solu- 
tion. Moreover, calculating the error on the solution of the linearized system 
does not provide information from a numerical point of view (sections 4.3 and 
4.6). Our purpose is not to verify the second order accuracy of the numeri- 
cal method (it has already been verified for the Shallow Water applications in 
[8, 4, 22]). But these errors computations allow us to verify that for these config- 
urations the nonlinear system has a linear behavior. And this linear behavior is 
perfectly captured by the scheme and the physical model. Precision is increased 
by increasing the order of the numerical method (as shown on Figure 6b and 
as illustrated in tables 1 and 2). The theoretical order of the schemes is not 
recovered for the wave test in section 4.3 (0.78 instead of 1 and 1.15 instead of 
2) this is due to the regularity of the solution of the initial system. Indeed, this 
solution is only continuous. In addition, the system has no term of order 2, then 
the interest in using a second order method is to save computing time for a given 
accuracy (see Figure 15). The semi-implicit treatment and the apparent topog- 
raphy give closed results. However, when friction increases (C/ = 0.005053), 
the semi-implicit method is closer to the linearised solution. 
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Figure 15: The wave "damping": The efficiency curves for C/ = 0.000202. 
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tcpu [s] 


\\Q - Qea; Ll 


tcpu [s] 


\\Q - QexWL^ 


tcpu [s] 


50 


0.314E-7 


0.49 


0.104E-7 


2.63 


O.llE-7 


2.7 


100 


0.159E-7 


1.93 


0.508E-8 


10.49 


0.546E-8 


10.7 


200 


0.801E-8 


7.7 


0.238E-8 


40.93 


0.263E-8 


42.6 


400 


0.397E-8 


30.93 


0.103E-8 


163.47 


0.118E-8 


170.32 


800 


0.192E-8 


123.28 


0.36E-9 


654.64 


0.433E-9 


681.26 


Regression 


y=-1.007x-13.32 




y=-1.2x-13.59 




y=-1.155x-13.72 





Table 2: The wave "damping": errors on the discharge Q and CPU times 
tcpu for Cf = 0. 

5 Conclusion 

In this paper we have considered the classical 1-D model of flow in an artery, 
we have presented a new numerical scheme and some numerical tests. The 
new proposed scheme has been obtained in following recent advances in the 
shallow- water Saint- Venant community. This community has been confronted to 
spurious flows induced by the change of topography in not well written schemes. 
In blood-flows, it corresponds to the treatment of the terms due to a variating 
initial shape of the artery, which arises is stenosis, aneurisms, or taper. To write 
the method, we have insisted on the fact that the effective conservative variables 
are: the artery area A and the discharge Q, which was not clearly observed up 
to now. The obtained conservative system has been then discretized in order to 
use the property of equilibrium of "the man at eternal rest" analogous to the 
"lake at rest" in Saint- Venant equations. If this property is not preserved in 
the numerical stencil, spurious currents arise in the case of a variating vessel. 
For sake of illustration, if the terms are treated in a too simple way, we have 
exhibited a case with an aneurism which induces a non zero flux of blood Q. 
In a pulsatile case, the same configuration creates extra waves as well. An 
analogous of dam break (here a tourniquet) was used to validate the other terms 
of the dicretisation. Other less demanding examples have been performed: linear 
waves without and with damping in straight tubes. Good behaviour is obtained 
in those pertinent test cases that explore all the parts of the equations. 

The well-balanced finite volume scheme that we propose preserves the vol- 
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ume of blood and avoids non physical case behavior. Thus we get a stable 
scheme and the accuracy is improved thanks to a second order reconstruction. 
The numerical method has been tested on examples which are more " gedanken" 
than physiological. This means that we only did a first necessary step. The next 
one is now to test more complex cases such as bifurcations, special boundary 
conditions and so on, in order to confront and fit to real practical clinical cases 
one of them being the case of aneurisms, in which as well we can introduce 
a source term in the mass equation in the case of leaks, we conclude that to 
look at those difBcult problems, a careful treatment of the numerical scheme is 
important. 

We thank ENDOCOM ANR for financial support, we thank Marie-Odile 
Bristeau and Christophe Berthon for fruitful discussions. 
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